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THREE-DIMENSIONAL ELASTIC STRESS AND DISPLACEMENT ANALYSIS 
OF FINITE CIRCULAR GEOMETRY SOLIDS CONTAINING CRACKS 
by John P. Gyekenyesi, Alexander Mendelson, and Jon Kring 
Lewis Research Center 

SUMMARY 

Displacement and stress distributions are calculated in finite circular bars, each 
containing a penny- shaped crack and loaded normal to the crack surface. Similar re- 
sults are presented for an annular plate containing internal, traction-free surface 
cracks. The Navier- Cauchy equations of elastic equilibrium are reduced to three sets 
of coupled, simultaneous, ordinary differential equations whose solutions are obtained 
along sets of lines in a discretized region. Since decoupling of these equations and their 
boundary conditions is not possible, a successive approximation procedure is applied to 
obtain their analytical solution. These analytical solutions are compared with known so- 
lutions, and the results of the circular bar are used to examine the rate of convergence 
and the accuracy of this method. The results obtained show that the method of lines 
presents a systematic approach to the solution of some three-dimensional elasticity 
problems with arbitrary boundary conditions. The advantage of this method over other 
numerical solutions is that good results are obtained even from the use of a relatively 
coarse grid. 


INTRODUCTION 

The object of a problem in elasticity is usually to calculate the displacement and 
stress distributions in an elastic body which is subject to given body forces or surface 
conditions. These distributions are the solutions of the applicable field equations which 
mathematically model the behavior of engineering materials. The solution of this gen- 
eral system of equations is, however, usually too difficult to evaluate. At present, few 
analytical solutions of three-dimensional problems exist (ref. 1), and even these solu- 
tions frequently require symmetry conditions to simplify the governing equations. Re- ' 
cently, with the introduction of large digital computers, the use of a number of approx- 


imate methods was attempted, but these methods yielded only partial results for these 
problems. Among these methods are the finite difference (ref. 2), the direct potential 
(ref. 3), the finite element (ref. 4), the eigenfunction expansion (ref. 5), and the line 
method of analysis (ref. 6). Of all these solution techniques, the line method of analysis 
is probably the least known and least used method in three-dimensional elasticity. Al- 
though the concept of this method for solving partial differential equations is not new 
’(ref. 7), its useful application in the past has been limited to simple examples (ref. 8). 

The line method lies midway between completely analytical and discrete methods. 

• The basis of this solution technique is the substitution of finite differences for the deriva- 
tives with respect to all the independent variables except one, with respect to which the 
derivatives are retained. This approach replaces a given partial differential equation 
with a system of simultaneous ordinary differential equations whose solutions can then 
be obtained in closed form. These solutions describe the dependent variable along lines 
which are parallel to the coordinate in whose direction the derivatives were retained. 
Because of their practical importance and inherent singularities, the work in this report 
is concentrated on the elastic analysis of three-dimensional, finite geometry solids 
which contain traction-free flaws or cracks. If we assume that the method of lines can 
be successfully applied to these solids with inherently large stress and strain gradients, 
its use for less complicated elasticity problems should present little difficulty. 

Early three-dimensional solutions of crack problems usually described the stresses 
near circular or ellipsoidal cavities enclosed in infinitely large solids. In reference 9, 
Sneddon applied Hankel transform methods to Love’s biharmonic strain function and re- 
duced the mixed boundary value problem of the axisymmetric half space to a set of dual 
integral equations. Only limited work has been done on three-dimensional solutions of 
crack problems in finite geometry solids . After making certain assumptions on the 
nature of the thickness variation of the stresses, Hartranft and Sih were successful in 
obtaining partial results for cracks in finite thickness bodies using variational methods 
(ref. 10) and eigenfunction expansions (ref. 5). In 1970, Cruse and VanBuren (ref. 11) 
used the direct potential method to analyze the stresses and displacements in a fracture 
specimen with a single edge crack. 

This report presents a simple and systematic approach to the elastic analysis of 
three-dimensional finite and circular geometry solids which contain traction-free cracks. 
A discussion of the elastic solutions, obtained by the same method of lines, of similar 
rectangular geometry solids can be found in references 6 and 12. The advantage of this 
method over other numerical techniques is that good results are obtained by using rela- 
tively coarse grids. This use of a coarse grid is permissible because parts of the solu- 
tions are obtained in terms of continuous functions. Additional accuracy in normal 
stress distributions is derived from the fact that they are expressed as first-order deri- 
vatives of the displacements and these derivatives can be analytically evaluated. Inher- 
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herently inaccurate numerical differentiation is required only for the evaluation of the 
shear stresses, but this presents no important loss of accuracy since they are an order 
of magnitude smaller than the normal stresses. For problems with geometric singular- 
ities, additional accuracy is derived from using a displacement formulation since the 
resulting deformations are not singular. 


SYMBOLS 


A(Kp» P) 

A ii 

a 

B i 

b 


c 

d 

E 


F, 


coefficient matrix of first-order differential equations, p = r, 9, z 
partitioned submatrices of e A ^, i = j =1,2 and p = 6, z 
crack radius 

particular integrals of first-order differential equations, i = 1, 2, . . . , 6 
outside radius 

surface crack length emanating from internal hole 
depth of surface crack along internal hole surface 
Young's modulus of elasticity 
dilatation or natural logarithm base 

initial value vectors for the first-order radial differential equations, 
i = 1,2 


G 

h r’V h z 

I 

K I 

Kr,K 0 ,K z 

L 

l 

m 

NR, N0,NZ 

n 

R 

r 0 


shear modulus of elasticity 

increments along cylindrical coordinate axes 

identity matrix 

stress intensity factor for opening mode 

coefficient matrices of second-order differential equations 

half-length of cylinder 

number of lines in radial direction 

number of lines in circumferential direction 

number of lines in given plane 

number of lines in axial direction 

distance from crack edge 

internal hole radius 
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r,e,z 

r(r) 

b(6) 

t 

to z) 

u 

v' 

w 

T? 

0 O 

X 

V 


cylindrical coordinates 

coupling vector for radial second-order differential equations 
coupling vector for circumferential second-order differential equations 
half-plate thickness 

coupling vector for axial second- order differential equations 

radial displacement 

circumferential displacement 

axial displacement 

variable of integration 

circumferential length of annular plate cutout 
Lame’s constant 
Poisson's ratio 


a r’ a 6’ °z 1 

„ „ r components of the stress tensor in cylindrical coordinates 

CT r0’ CT rz’ °Qz) 


°U 


matrizant of A(K r ,r) 

partitioned submatrices of the matrizant of A(K r , r) 
Laplacian 


The symbols used in the appendixes are defined when they are introduced. 


REDUCTION OF NAVIER -CAUCHY EQUATIONS TO SYSTEMS 
OF ORDINARY DIFFERENTIAL EQUATIONS 

The fifteen individual field equations of linearized elasticity may be combined to 
form the three Navier- Cauchy equations of elastic equilibrium which in cylindrical coor- 
dinates are written as 


— + (1 - 2v) 
dr 




u - 


2 dv 
2 dd 


= 0 


( 1 ) 


4 



li§+ (1 . 2v) 
r 3 0 



v + 


2 Su 

2 30 


= 0 


( 2 ) 


— + (1 - 2i/)V 2 w = 0 
3z 


(3) 


where the body forces are assumed to be zero and the dilatation is 


„ 3u 1 3v , u , 3w 
e - — + + — + — 

3r r 30 r 3z 


(4) 


The stress-displacement relations, obtained by substituting the strain-displacement 
relations into Hooke's law, can be expressed in the following form: 


CT r = Xe + 2G 


3u 

3r 


(5) 


Uq — Xe + 2G 


/ 1 3v | u 
\r 30 r, 


( 6 ) 


cr = Xe + 2G — 
z 3z 


(7) 


_ r / 1 3u v 3v 

r0 - " _ + ~ 

\r 30 r Sr, 


( 8 ) 


CT rz = G + — 
rz \3r 3z 


(9) 


_ r /3v , 1 3w 


( 10 ) 


For a finite geometry body with circular boundaries, we construct three sets of 
parallel lines in the direction of the coordinates as shown in figure 1. Approximate 
solutions of equations (1) to (3) can then be obtained by developing solutions of ordinary ' 
differential equations along the radial, circumferential, and axial lines, respectively. 
For equation (1), we define the displacements along the radial lines as u^Ug, . . 

The derivatives of the circumferential displacements on these lines with respect to 0 


5 



Figure 1. - Sets ot lines in direction of cylindrical coordinates. 


are defined as v| • • . , v| and the derivatives of the axial displacements with 

respect to z are defined as w|p w| g, • • . , w|^. These displacements and deriva- 
tives can then be regarded as functions of the radius only, since they are variables on 
radial lines. If these definitions are used, the ordinary differential equation along a 
generic, radial line ij (a double subscript is used here for simplicity of writing) of fig- 
ure 1 may be written as 


d^u.. du-- u.. 
1] ! 1 1] _ _ii 


dr 2 r dr r 2 


(i - M 
2(1 - V ) 


— + — ) u ii + — ( u i+ l i + u i 1 i) 

iA* d 11 A * ,+1 ' 1 Il>J 




f ii (r) 
2(1 - v) 


= 0 (ID 
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where 


*1 2 ^ r dr 

r 


i] 


dw| 

dr 


i] 


( 12 ) 


and 


v 


dv 

de 



Similar differential equations are obtained along the other radial lines. Since each equa- 
tion has the terms of the displacements on the surrounding lines, these equations consti- 
tute a system of ordinary differential equations for the displacements UpUg, . . . , u^. 

Noting that a second-order differential equation can satisfy only a total of two 
boundary conditions and since three-dimensional elasticity problems have three boundary 
conditions at every point of the bounding surface, some of the boundary data must be in- 
corporated into the surface line differential equations. Hence, conditions of normal 
stress and displacement are enforced through the constants of the homogeneous solutions 
while shear stress boundary data must be incorporated into the differential equations of 
the surface lines. It is the application of the specified shear conditions that permits the 
use of central difference approximations when surface line differential equations are 
constructed. For the first radial line of figure 1, the use of zero shear stress boundary 
conditions in the radial direction on the r, z and r, 0 coordinate planes gives, respec- 
tively, the following imaginary radial line displacements: 


= u 0 + 2h„r — 


I0e " u 2 + au e 


dr 


2 V 1 


u 10n _U N0+1 + 2h z^T 


(13) 


Equations (13) must then be used in the application of central difference approximations 
when the ordinary differential equation for the first radial line is generated. Additional 
details on the construction of these equations can be found in reference 6. It is conven- 
ient to nondimensionalize these equations with respect to some characteristic dimension. ' 
Hence, we introduce the following variables: 
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~ u 
u = - 

a 

T = — 

a 

>| * 
ii 



~ y 

V = - 

a 

0=0 

**0 _ h 0 

► (14) 

* 

3? 

ii 

13 

'z = - 
a 

N 

II 

» 1*°" 
v 


where a is the 
the radial lines 

crack radius . If matrix notation is used, 
can be expressed as 

the differential equations along 



1 d 
r dr 



l x 1 


rsj rsj r*j 

= K r (r) u + r(r ) 
unxi i x i 


(15) 


where the coefficient matrix K r (r ) is a function of the radius, the coordinate incre- 
ments h 0 and h z , and Poisson’s ratio. 

In a similar manner, for the solution of equations (2) and (3) ordinary differential 
equations are constructed along the circumferential and axial lines, respectively. 
These equations, however, are different in form from equation (15) and may be written 
as 


^ = K. V s (?) 
d? 2 

m x l mxmmxlmxl 


, 2 ~ 
d w 


dz 
n x 1 


K z w + t(z) 
nXnnxl nx 1 


(16) 


(17) 


These two sets of equations are linear, second-order differential equations with con- 
stant coefficients. The coupling terms in each set of differential equations are grouped 
into the vectors r(r ), s(o ), and t(z: ), which are the nonhomogeneous terms in the pre- 
vious equations. Since the sum of the elements in any given row of the K r (r ) and K z 
coefficient matrices is zero, they are both singular. In addition, we find that is 
also singular, although this is not as evident from its elements as in the case of the other 
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coefficient matrices. A detailed listing of equations (15) to (17) may be found in refer- 
ence 6. 

SOLUTION OF THE SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 

The details of a solution technique for the type of simultaneous ordinary differential 
equations given by equations (15) are discussed in appendix A. For the case r * 0, the 
solution vector for the radial displacement is 

u(r) = 4 « n (r) \F 1 + Bj(r)]+i fi 12 (r) [f 2 + B 2 (r)] ( 18 ) 

Z x 1 l x l l x 1 Z x 1 lx l lx 1 Z x 1 

where 

Fpru 

r= initial 


r r initial 

and fij^(r ) and j 2 ( r ) are infinite matrix integral series whose definitions and evalu- 
ation are discussed in appendix A. Vectors B^(r)and B 2 (r ) represent particular 
solutions of equations (15), and they are given by 



(V ^ 

j_(r ) = - ^11^ 12^22 " ^21^11^12) *^77 


' initial 





- ^ 21 ^ 11^12 



( 20 ) 


while vectors and F 2 are initial value vectors whose evaluation from given bound- 
ary conditions is discussed in appendix C. Differentiating equation (18) with respect to 
r yields 
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( 21 ) 


/"o n*j 


u(r) = 


1 


r2i^ r) 

r L 

Z x 1 Z x Z ZxZZxl Z x 1 


[ F i +B i< 7 )] 


) 7 ^I2( r ) 


22 
Z x Z 


r 2 


[f 2 + B 2 (r)] 


Z X Z zxl Z X 1 


where fi 24 (r ) an< ^ ^22^ r ) are ma ^ r i x integral series similar to ^^(r ) and £Zj 2 (r )• 
For differential equations of the type given by equations (16) and (17), a detailed and 
analogous solution procedure is discussed in appendix B. As a result of this develop- 
ment, the circumferential displacement vector may be written as 

v(e ) = A 11 (e)[v(o) + b 3 (0)] + a 12 (0)[v(o) + B 4 (e)] (22) 

where A 4 j(6i ) and A^ 2 ( 0 ) are matrix functions: 

A n (0) = cosh kJ/ 2 0 (23) 

A 12 (0 ) = YL^ 2 sinh K 1 / 2 0 (24) 

The matrix K^ 2 is defined as a matrix whose square is equal to K 0 . Vectors v(0) 
and v (0) are initial value vectors whose evaluation from given boundary conditions is 
discussed in appendix C. Vectors Bg(^)and B^( 0 ) represent particular solutions of 
equations (16) and are given by 


B 3 (0) 




A 11 (tj)s(tj) dr] 


Differentiating equation (22) with respect to 0 gives 

v (' 0 ) = K 0 A 12 ( 0 )[v (0) + B 3 ( 0 )] + A n ( 0 )[v (0) + B 4 ( 0 )] 


(25) 


(26) 


( 27 ) 


A solution vector similar to equation (22) may be constructed for equations (17) which 
defines the axial displacements. Since r(r)) in equations (19) and (20) is unknown, we 
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start the solution of the problem by assuming values of v, dv/dr , dw/dr, v, dv/dr, 
and dw/dr along the radial lines. The initial values of these required quantities were 
taken to be zero in this report. Then equations (18) and (21) will give the first estimate 

rv /I \ tv { 1 \ 

for the vectors u(r) v ' and u(r) v It is assumed, of course, that vectors and 
F 9 are known or can be found from given boundary conditions. Using these calculated 

ro rsj ( 1 \ /v /v \ f 1 ^ rsJ . ( 1 } 

values of u (r ) v ' and u (r y \ we can evaluate the vector s ( 6 p ' where we must use 
similarly assumed values of dw/d6> and dw /dd along the circumferential lines. 


Equations (22). and (27) give the first values of v ( 6 )^ and v (0 ) K±> . First values of 
w ( z ) v ’ and w ( z p ’ can then be calculated by using the first estimates of the radial and 
circumferential displacements and their derivatives in the coupling vector t( z )^. At 


r*j r>j 


(1) 


this point we return to equations (18) and (21) and calculate the second values of u (r ) 

As (O } 

and u (r ) v ' based on the first values of the circumferential and axial solutions 


( 2 ) 


If the values of u(r)^, u(r)^, v(6>)^, v(0)^, w(z)^, and w (z )^ converge 
with the repetition of this procedure, an approximate solution of a given problem will be 
determined. In general, the convergence rate for these successive calculations is de- 
pendent on the accuracy to which the matrix functions A.j and the matrix integral ser- 
ies can be evaluated. Sufficiently large errors in these matrix variables will lead 
to divergence in the successive approximation calculations. 

Since the coupling vectors r(r ), s(6>)> and t(*z ) involve displacements and their 
derivatives defined only at the nodes, finite difference calculus must be used in evaluat- 
ing their elements. Hence, all the particular integrals are calculated by a suitable 
numerical integration technique. 

Once the displacement field in the bar has been calculated and the successive ap- 
proximation procedure has converged, the normal stress distributions along the three 
sets of parallel lines can be obtained from 


a r _ (X + 2G)(u) Along rad . al + X 


lines 


fi-L* 

\ r r 


/v 

v + w 




/ 


Along radial 
lines 


(28) 


o e = (X + 2G) 


(B) 


r*j , 


+ X(u + w)^i on g c , rcum f er . 

Along circumfer- ential Unes 

ential lines 


(29) 


o 2 = (x + 2G)(W ) Along 
lines 


Along axial 
lines 


(30) 
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Note that these equations involve only derivatives that can be analytically evaluated. The 
shear stresses at each node can be calculated from equations (8) to (10), but finite dif- 
ference approximations must be used for the required displacement gradients. 


NUMERICAL RESULTS 

Solid Cylindrical Bar with Penny-Shaped Crack 

Figure 2 shows a cylindrical bar containing a penny-shaped crack and loaded by a 
uniform normal stress distribution. For problems with axisymmetric geometry, the 
circumferential displacement is inherently zero at every point and all the remaining 
variables are independent of the circumferential coordinate 6 . The two sets of parallel 
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lines needed for the solution of this problem are also shown in this figure. Note that the 
crack edge is assumed to be midway between adjacent nodes, specifying normal stress 
and displacement boundary conditions, respectively. 

The solution of this problem was obtained by using two different sets of lines along 
the coordinate axes so that the convergence of the finite difference approximations could 
be checked. The successive approximation procedure required for decoupling the two 
sets of ordinary differential equations was terminated when the difference between sue- 
cessively calculated displacements at every point was less than a preset value (10 ). 

Selected results from these calculations are shown in tables I to V and figures 3 to 7. - 
For easy comparison of data, some of these figures include Sneddon’s results for an 


TABLE I. - NONDIMENSIONALIZED RADIAL DISPLACEMENTS Eu/a Q a 
FOR SOLID CYLINDRICAL BAR WITH PENNY-SHAPED CRACK 
UNDER UNIFORM NORMAL TENSION 


(a = 1. 0, b = 1. 77, L = 1. 68 (16 axial and 13 radial lines)!] 


'z 

r 


0.000 

0.235 

0.471 

0.706 

0.941 

1.06 

1.294 

1.530 

1.770 

0.00 

0.000 

-0. 168 

-0.333 

-0.488 

-0.579 

-0. 609 

-0. 653 

-0.705 

-0.773 

.28 



-.080 

-. 153 

-.217 

-.297 

-.371 

-. 529 

-.639 

-.725 

.56 



-.045 

-.091 

-. 149 

-.239 

-.300 

-.433 

-.549 

-.646 

.84 



-.038 

-.083 

-. 142 

-.223 

-. 273 

-.379 

-.482 

-.580 

1. 12 



-.039 

-.085 

-.142 

-.214 

-.255 

-.344 

-.435 

-.530 

1.40 



-.034 

-.074 

-. 124 

-. 188 

-.224 

-.303 

-.388 

-.483 

1.68 



-.006 

-.023 

-.058 

-. 114 

-. 150 

-.233 

-.325 

-.425 


infinite solid (ref. 9). The data of figure 3 clearly show the advantage of the line method 
over other numerical solutions. A relatively coarse grid of nine axial and nine radial 
lines gave almost identical results to those obtained by using a 16 by 13 grid. Since the 
bar is of finite size, the crack opening displacement is expected to be higher than 
Sneddon's solution. Consistency of the results with this conclusion is obvious from 
figure 6. It is noteworthy that the results correspond to elliptical crack profiles in all 
cases. The maximum dimensionless crack opening is plotted for several crack to cyl- 
inder radius ratios in figure 7. The data from figure 6 match the results of Sneddon 
and Welch (ref. 13) very closely, and, as expected, the maximum crack opening is 
slightly higher for finite length cylinders. 

Figure 5 shows the stress distribution normal to the crack plane as a function of the 
distance from the crack edge . As shown by Sneddon for an infinite solid, this stress 
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TABLE II. - NONDIMENSIONALIZED AXIAL DISPLACEMENTS 


Ew/ff Q a FOR SOLID CYLINDRICAL BAR WITH 
PENNY-SHAPED CRACK UNDER 
UNIFORM NORMAL TENSION 


[a = 1.0, b = 1. 77, L = 1. 68 (16 axial and 13 radial lines)!] 


r 

z 


0.000 

0.280 

0.560 

0. 840 

1. 120 

1. 400 

1.680 

0.000 

1.394 

1.544 

1.646 

1.766 

1.928 

2. 126 

2.329 

. 235 

1.351 

1.495 

1.595 

1.722 

1.892 

2. 097 

2.306 

.471 

1.219 

1.359 

1. 472 

1.622 

1. 812 

2.029 

2.250 

. 706 

.973 

1. 113 

1.273 

1.473 

1.696 

1.933 

2. 167 

.941 

.515 

.746 

1.029 

1.302 

1.565 

1.823 

2.072 

1. 060 

.000 

.542 

.913 

1.221 

1. 502 

1. 768 

2.024 

1.294 



.384 

.758 

1.090 

1.390 

1.670 

1.937 

1.530 



.334 

. 673 

.999 

1.304 

1.592 

1. 866 

1.770 



.291 

.610 

.935 

1.245 

1.535 

1. 812 


TABLE IH. - NONDIMENSIONALIZED RADIAL STRESS a r /° 0 FOR 
SOLID CYLINDRICAL BAR WITH PENNY-SHAPED CRACK 
UNDER UNIFORM NORMAL TENSION 


[a =1.0, b = 1. 77, L =1. 68 (16 axial and 13 radial lines)!] 


z 

7 


0. 000 

0.235 

0.471 

0.706 

0.941 

1.060 

1.244 

1.530 

1.770 

0.00 

-1.060 

-1.043 

-1.017 

-0.883 

-0.522 

1. 134 

0.376 

0.116 

0.000 

.28 

-.479 

-.451 

-.375 

-.270 

-.262 

-. 138 

-.095 

-.017 



.56 

-. 148 

-. 141 

-. 123 

-. 124 

-. 155 

-. 150 

-. 113 

-.044 



.84 

.0211 

.0144 

.002 

-.0251 

-.058 

-.065 

-.056 

-.026 



1. 12 

. 127 

. 116 

.093 

. 058 

.023 

. 009 

-.004 

-.003 



1.40 

.246 

.230 

. 193 

. 142 

.088 

. 064 

.028 

.012 



1.68 

.470 

. 429 

.352 

.247 

. 136 

. 085 

.008 

-.021 







TABLE IV. - NONDIM E NSIONA LI Z E D CIRCUMFERENTIAL STRESS a 0 /cr o 
FOR SOLID CYLINDRICAL BAR WITH PENNY-SHAPED CRACK 
UNDER UNIFORM NORMAL TENSION 


[a = 1. 0, b = 1. 77, L = 1. 68 (16 axial and 13 radial lines)]] 


z 

r 


0.000 

0. 235 

0. 471 

0.706 

0.941 

1. 060 

1. 294 

1.530 

1.770 

0. 00 

-1.060 

-1. 061 

-1.047 

-0.986 

-0. 790 

0.909 

0. 125 

-0. 020 

-0. 108 

.28 

-.479 

-.461 

-.401 

-. 291 

-.112 

. 108 

.015 

-.023 

-.051 

.56 

-. 148 

-. 138 

-. 106 

-.055 

. Oil 

. 049 

. 037 

. 022 

.024 

.84 

.021 

.023 

.032 

. 047 

.062 

. 067 

.062 

. 055 

. 057 

1. 12 

. 127 

. 124 

. 120 

. 113 

. 104 

. 099 

. 086 

. 075 

.064 

1.40 

.246 

.23 

.218 

. 191 

. 160 

. 144 

. 116 

. 096 

.072 

1.68 

. 470 

. 45 

.401 

.333 

.257 

. 220 

. 154 

. 112 

.089 


TABLE V. - NONDIMENSIONALIZED AXIAL STRESS u z / o q FOR 
SOLID CYLINDRICAL BAR WITH PENNY-SHAPED CRACK 
UNDER UNIFORM NORMAL TENSION 


[a = 1. 0, b = 1. 77, L = 1. 68 (16 axial and 13 radial lines)]] 


z 

r 



0.000 

0.235 

0. 471 

0 

706 

0 

941 

1 

060 

1 

294 

1. 

530 

1. 770 

0 

00 

0.000 

0.000 

0. 000 

0 

000 

0 

000 

3 

320 

1 

512 

1. 

206 

0.989 


28 

.095 

.093 

. 146 


319 


874 

1 

513 

1. 

365 

1. 

201 

1. 082 


56 

.277 

.295 

.386 


592 


950 

1 

150 

1 

230 

1. 

186 

1.172 


84 

.515 

.542 

.624 


770 


957 

1 

041 

1 

122 

1. 

136 

1. 159 

l 

12 

.736 

.759 

.808 


885 


972 

1 

010 

1 

059 

1. 

082 

1. 092 

l 

40 

.900 

.915 

.933 


960 


990 

1 

003 

1 

022 

1. 

036 

1. 037 

l 

68 

.999 

.998 

.999 


999 

1 

000 


999 


995 


993 

. 990 






grid). 



(b) Dimensionless axial displacement distribution <9 by 9 grid). 

Figure 3. - Dimensionless axial displacement distribution 
for solid cylindrical bar with penny-shaped crack. 


distribution should approach infinity near the crack edge as the inverse square root of 
the distance from the crack edge. Establishment of this type of singularity is, however, 
difficult when numerical methods are used because values of the normal stress are 
needed within a distance of 0. 05 a or less from the crack tip. For the range of r 
shown in figure 5, this inverse square root singularity is not defined. However, for the 
range shown, the obtained stress curve closely resembles Sneddon's solution (ref. 15). 
As expected, the absolute value of this stress is greater for a finite size bar than for 
Sneddon’s infinite solid. Tables I to V show selected results from the computer listings. 
The accuracy of the normal stress and displacement boundary conditions can easily be 
noted from the numerical data listed. 
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grid). 



(b) Dimensionless radial displacement distribution (9 by 9 grid). 

Figure 4. - Dimensionless radial displacement distribution 
for solid cylindrical bar with penny-shaped crack. 



Figure 5. - Dimensionless axial stress distribution (16 by 13 grid) 
for solid cylindrical bar with penny-shaped crack at z - 0. 
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Figure 6. - Dimensionless crack opening displacements for 
solid cylindrical bars with penny-shaped cracks of various 
lengths and radii. 



Figure 7. - Dimensionless maximum crack opening as 
function of crack to cylinder radius ratio. 


Annular Plate with Internal Surface Cracks 

As a first attempt at the solution of a general three-dimensional problem in cylin- 
drical coordinates, the problem of figure 8 was investigated. In order to minimize the 
numerical computations, we have assumed four symmetrically located internal surface 
cracks. Because of the symmetric geometry and loading, only one-sixteenth of the 
original plate has to be discretized. Figure 9 shows this region of interest and the as- 
sumed crack geometry. Note that nondimensionalization with respect to the outside 
radius is more convenient in this case since the crack has two characteristic dimensions. 

Displacement and normal stress distributions were calculated along the lines of 
figure 9 using a grid of 16 lines in all directions. Because of the relatively coarse grid 
involved, the results appear only in tabular form (tables VI to XI). It must be noted, 
however, that because of the unknown nature of the resulting solutions, the use of a 
coarse grid is always recommended in generating the first set of displacements. Since 
the construction of a general computer program for this problem requires a great 
amount of effort, the results in tables VI to XI were calculated by using only one set of 
lines with NR = N0 = NZ = 4. 

Tables VI to VIII contain the dimensionless displacement distributions. These re- 
sults show that below the crack plane the circumferential displacements are essentially 
zero while the maximum crack opening is at 9 = z =0 and r = 0.25. Although both 
radial and circumferential displacement fields are extensional, the axial displacements 
are negative, which indicates a contraction in that direction. The calculated normal 
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TABLE VI. - DIMENSIONLESS RADIAL DISPLACEMENTS 


Eu/CT Q b FOR ANNULAR PLATE WITH INTERNAL 
SURFACE CRACKS UNDER UNIFORM RADIAL 


TENSION ON OUTSIDE SURFACE 


0, 

deg 

r 

0.25 

0.50 

0.75 

1.00 

z = 0. 00 

0 

0.435 

0.671 

0.786 

0.930 

15 

.940 

. 890 

.896 

.998 

30 

.168 

.992 

.968 

1.058 

45 

1.240 

1.026 

.994 

1.081 

z = 0. 10 

0 

0.506 

0. 680 

0.789 

0.930 

15 

.907 

.860 

.886 

.994 

30 

1.092 

.953 

.957 

1.055 

45 



1. 152 

.986 

.984 

1.076 

z =0.20 

0 

0.724 

0.656 

0.779 

0.927 

15 

. 795 

.761 

.856 

.965 

30 

. 898 

. 856 

.932 

1.048 

45 

.942 

. 892 

.961 

1.073 

z =0.30 

0 

0.716 

0. 659 

0.769 

0.915 

15 

.752 

.727 

. 845 

.978 

30 

. 815 

. 809 

.926 

1.049 

45 

. 845 

. 843 

.957 

1.077 




TABLE VII. - DIMENSIONLESS CIRCUMFERENTIAL 


DISPLACEMENTS Ev/a Q b FOR ANNULAR PLATE 
WITH INTERNAL SURFACE CRACKS UNDER 
UNIFORM RADIAL TENSION ON 


OUTSIDE SURFACE 


z 

0, deg 

0 

15 

30 

45 

r =0.25 

0.00 

0.718 

0.566 

0.302 

0. 000 

. 10 

.592 

.450 

.238 



.20 

.000 

.073 

.063 



.30 

. 000 

.009 

.012 



r =0.50 

0.00 

0.425 

0.285 

0. 138 

0. 000 

. 10 

.350 

.232 

.115 



. 20 

. 000 

.075 

.057 



CO 

o 

. 000 

.023 



.026 





r =0.75 

0. 00 

0. 000 

0.037 

0.029 

0»000 

. 10 



.033 

.026 



.20 



.023 

.020 



.30 



.014 

.014 



r =1.00 

0.00 

0. 000 

0. 006 

0. 004 

0. 000 

. 10 



.007 

.005 



. 20 



.009 

.007 



.30 



.010 

.007 






TABLE Vm. - DIMENSIONLESS AXIAL DISPLACEMENTS 


Ew/a Q b FOR ANNULAR PLATE WITH INTERNAL 
SURFACE CRACKS UNDER UNIFORM RADIAL 


TENSION ON OUTSIDE SURFACE 


r 

z 


0.00 

0 . 10 

0.20 

0.30 



0 = 0 ° 

0.25 

0 . 000 

- 0.218 

- 0.342 

- 0.416 

.50 



-. 108 

-.217 

-.322 

.75 



-.089 

-. 175 

-.256 

1.00 



-.076 

-. 152 

-.226 

0 = 15 ° 

0.25 

0 . 000 

- 0 . 098 

- 0.216 

- 0.336 

.50 



-.049 

-. 131 

-.232 

.75 



-.075 

-. 154 

-.237 

1.00 



-.074 

-.149 

-.226 

0 = 30 ° 

0.25 

0.000 

- 0.068 

- 0. 167 

- 0.281 

.50 



-.040 

-.106 

-. 195 

.75 



-.067 

-. 141 

-.222 

1.00 



-.075 

-. 152 

-.231 

0 = 45 ° 

0.25 

0.000 

- 0 . 061 

- 0. 153 

- 0.262 

.50 



-.039 

-. 101 

-.186 

.75 



-.065 

-. 137 

-.218 

1.00 



-.075 

-. 153 

-.233 




TABLE IX. - DIMENSIONLESS RADIAL STRESS 


DISTRIBUTION <r^/ <7 Q FOR ANNULAR 
PLATE WITH INTERNAL SURFACE 
CRACKS UNDER UNIFORM 
RADIAL TENSION ON 


OUTSIDE SURFACE 


0, 

deg 

r 

v 0.25 

0.50 

0.75 

1. 

00 

z = 0. 00 

0 

0.000 

0.020 

1.062 

1.000 

15 



.346 

. 837 



30 



.293 

.748 



45 



.282 

.726 



z =0.20 

0 

0.000 

0.027 

1.048 

1.000 

15 



. 405 

. 872 



30 



. 424 

.801 



45 



.424 

.783 


I 

z = 0. 10 

0 

0.000 

1.305 

1. 026 

1.000 

15 



.968 

.980 



30 



. 840 

.932 



45 



. 801 

.915 



z =0.30 

0 

0.000 

0. 872 

1.005 

1.000 

15 



.950 

1.004 



30 



.966 

.986 



45 



.956 

. 975 






TABLE X. - DIMENSIONLESS CIRCUMFERENTIAL 


STRESS a„/a FOR ANNULAR PLATE WITH 

v O 

INTERNAL SURFACE CRACKS UNDER 
UNIFORM RADIAL TENSION ON 
OUTSIDE SURFACE 


0, 

r 

deg 






0.25 

0.50 

0.75 

1.00 

z = 0. 00 

0 

0.000 

0.000 

1.702 

1.251 

15 

.308 

.803 

1.542 

1.339 

30 

.209 

1.043 

1.460 

1.416 

45 

.240 

1. 145 

1.441 

1.443 

z = 0. 10 

0 

0.000 

0.000 

1.659 

1.246 

15 

.703 

.957 

1.534 

1.331 

30 

. 803 

1. 196 

1.472 

1.410 

45 

. 863 

1.270 

1.456 

1.439 

z = 0. 20 

0 

4.569 

2. 821 

1.506 

1.236 

15 

3.609 

2.011 

1.526 

1.313 

30 

2.971 

1.706 

1.506 

1.400 

45 

2.748 

1.629 



1.498 



1.432 

z =0.30 

0 

3.050 

1. 858 

1.442 

1.220 

15 

3. 110 

1. 879 

1.500 

1.301 

30 

3. 171 

1. 825 

1.523 

1.397 

45 

3. 167 

1.775 

1.524 

1.432 
































stresses are summarized in tables IX to XI. It is expected that the circumferential 
stress be the maximum tensile stress in the body since the crack plane is normal to this 
coordinate. The interaction between the hole and crack stress fields should also be most 
pronounced at the inside radius; from this we conclude that the stress should be maxi- 
mum at that point. The results of table X seem to confirm these observations, although 
we must also remember that the nodes are closer to the crack boundary below the crack 
surface than along the radius. 

Because of the use of a coarse grid, the numerical results for this example are 
somewhat inaccurate in magnitude but they do indicate some previously unknown varia- 
tions in the stress field for this problem. This conclusion is possible in that the line 
method does not usually require a fine grid for good results as was shown in the previ- 
ous section. These results also demonstrate that the method of lines permits the com- 
putation of the displacement and stress fields for a general three-dimensional problem. 


STRESS INTENSITY FACTOR 

It is customary in fracture mechanics to describe the plane elasticity crack opening 
displacement as a superposition of three basic deformation modes (ref. 14). In terms 
of the stress intensity factor for the opening mode Kj, the plane elasticity crack dis- 
placements near the crack tip are given by (ref. 14) 


ly=0 


_ 2 ( 1 -v) 



plane strain 


(31) 


1 y=° " (i + „)g 


/— plane stress 
2n 


(32) 


Since three-dimensional problems are neither in a state of plane strain nor plane stress, 
the definition of a stress intensity factor for these problems must be first established. 
The problem we consider in detail is Sneddon's penny-shaped crack solution. Refer- 
ence 15 gives the crack opening displacement as 


w 


z=0 


4cr Q (l - v ) 
jtE 


J ~ 2 2 

f a - r 


( 33 ) 


which for small values of R, where R = a - r, becomes 
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( 34 ) 


4cr 0 (l - if) , 

w|„_ n = — - V2aR 


1 z=0 


nE 


Note that adjustment of coordinates must be made between the equations quoted from 
references 14 and 15. Paris and Sih (ref. 14) list the stress intensity factor for this 
problem as 




(35) 


In terms of this stress intensity factor, the crack opening displacement (eq. (34)) 
becomes 


w z=0 = K I 


2(1 - 


E 


JL) Jl R _ 2(1 - v) Y JR 
In G 1 


2n 


(36) 


Rearranging this equation in terms of the known dimensionless displacements 


(E/a Q )(w/a) 


z=0 


gives 


E 


E w 
°0 a 


C I K I = 


z=0 


(37) 


where 


C - 4 (! - v 2 ) 
E V2?ra 


(38) 


Then a plot of equation (37) as VR/a — 0 gives the desired value of (E/<7 q)CjKj from 
which an equivalent stress intensity factor for finite geometry cylinders can be calcu- 
lated. Note from equation (36) that the definition of Kj given in equation (35) implies 
the plane strain crack opening displacement of equation (31). 

For problems in which the crack opening displacement varies in the thickness direc- 
tion, such as the annular plate with internal surface cracks, the stress intensity factor 
obtained previously will be a function of the thickness variable. However, if we were to 
account for the nonplane strain conditions near the surface by using equation (32) or a 
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corrected equation (31) for the definition of Kj, the stress intensity factor would become 
a constant across the thickness by definition. It should be noted that this description of 
Kj is completely arbitrary and its significance in three-dimensional elasticity theory is 
of doubtful value. However, values of Kj are still presented so that a comparison will 
be possible between the calculated results and the published plane strain solutions 
(ref. 14). 

Figure 10 contains the calculation of Kj according to equation (37) for the penny- 
shaped crack problems. The value of Kj calculated from Sneddon's solution agrees 



0 .2 .4 .6 .8 1.0 


Figure 10. - Calculation of the stress intensity factor Kj fora 
solid cylindrical ter with a penny shaped crack (16 by 13 grid). 

with equation (35), while for the finite size bar we obtain a value of Kj equal to 1. 485 
OQ-y/a. Hence, the finite bar discussed in figure 10 has an approximately 31 percent 
higher stress intensity factor than the infinite solid. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 22, 1973, 

501-21. 
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APPENDIX A 


SOLUTION OF SIMULTANEOUS DIFFERENTIAL EQUATIONS 
WITH VARIABLE COEFFICIENTS 


Noting that the differential operator of equations (15) can be written as d/dr[l/r x 
d/dr (ru)], the following variables are introduced in order to obtain a system of first- 
order differential equations: 


A 


Uj = r Up U 2 = r Uj, . . . , = r u . 


1 1 , 1 H ^ ~ . 1 H ~ . 

U *+l = ?dT rUl ’ U ^+2 = ~^ (ru 2 ) > * * • ’ U 2 1 = ~^ (rU Z } 


> (Al) 


In terms of these variables, equations (15) can be written as 


where 


dU 

dr 

21 x 1 


A(r) U r(r) 
21 x 21 21 x 1 21 x 1 


A(r) 

21 x 21 


0 

l x l 


^K r (r) 

r 

L l x l 


rl 

l x i 

0 

l x l m 


r(r) 
21 x 1 


' 0 ' 

l x 1 

r(r) 
l x 1 


(A2) 


(A3) 


(A4) 


Following the Peano- Baker form of integration (ref. 16), the solution of equation (A2) is 
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U(r) = nJ(A) U(0) + «J(A) ^[^(A)] ^(ij) d V 


(A5) 


21 x 1 21 x 21 21 x 1 2Z x 2Z 


21*21 21 x 1 


where vector U(0) consists of the boundary values of (r u) and [l/r d/dr (r u)3 at the 
initial value of r which is taken as zero in this case. The matrizant of A(r) is an 
infinite matrix integral series given by 

J2q(A)=I+ J A(Tj 1 )drj 1 + J A(t) 2 ) dr? 2 j ^ 2 A^j) drjj 


+ 


j £ A(tj 3 ) dr? 3 Jj 3 A(t} 2 ) dr7 2 J ^ 2 drjj + . . . (A6) 


Substituting matrix A(r ) into equation (A6) leads to the following four matrix integral 
series in terms of the coefficient matrix K(r ): 


n 11 (r)= I + 

lx l lx i 


r" 

f 

*/0 


t) 2 I d? ] 2 


r* ± 

Jo 


l x l 


Kj/tjj) di ) ] 


+ . 


I X l 


(A7) 


^12^ r ^ ~ ^ 
lx l lx l 


f 


Vl <ty] 


/ 

Jo 


f 


*3 j n 2 

1 ^3^3 / — K r^2^ d?7 2 / ^1 1 d,? l + 

. " 2 Jb ‘ 

Z x Z Z X l lx l 


(A8) 
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From the definitions of these matrix series we can easily conclude that their initial 
values are 
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W ll^ r )l r =0" 1 
n i2^^l r=0 = 0 

fi 21^lr=0 = 0 
fi 22^Hr=0 = I > 


(A 13) 


Since the matrix series (A7) to (A10) are usually difficult to evaluate by performing the 
indicated integrations, the previous simultaneous matrix differential equations may be 
solved by a suitable numerical technique instead. The necessary initial conditions for 
these numerical solutions are given by equations (A13). For the examples presented in 
this report, the single step Runge-Kutta integration method was used to solve equa- 
tions (All) and (A12). 

The particular integral of equation (A5) is written as 

^ 1 1 (f? ) ^12^ 

In evaluating the inverse of fi by partitioning, the partitioned integrals (A 14) are given 
by 



r 

B i(r) 

l x l 


> = 


Bo(r) 


B-. (r ) = 


-r 


^11^12^22 " ^21^11^12^ ^ (n ) d.77 


(A15) 


Bo(r ) 





)'■ 


r(r?) dr] 


(A16) 


If it is assumed that r(rj) is known, the previous particular integrals can easily be eval- 
uated. When the partitioned form of the matrices is used, equation (A5) may be written 
as 
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r u 


> = 




12 ' 


ft 


r ~ ~ ~\ 

r u 


f' r\ 
B, (r ) 




> +< 


(A17) 


Vs. 


1 d 

/VJ 

r dr 


(ru) 


^2 1( r ) fi 22^ r jJ 


1 1 /~~\ 
— (r u ) 

r dr 




r =0 


B 2 (r) 
L j 


7 


from which the solution vectors (18) and (21) can be directly constructed. 

It may be noted that for axisymmetric problems the following identities can be used 
to evaluate the particular integrals: 


v-1 


fi ll _ ( n 22 “ ^21^11^12) 

fi 12 =fi ll°12( fi 22 " ^21^11^12) J 


(A18) 


These identities are obtained from noting that the inverse of the matrizant Qq (A) is 
given by 



(A19) 


when the coefficient matrix is of the form (A3) and the matrix Kj. has constants for its 
elements. 
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APPENDIX B 


SOLUTION OF SIMULTANEOUS DIFFERENTIAL EQUATIONS 
WITH CONSTANT COEFFICIENTS 


For the solution of equations (16) we introduce the following new variables: 


V l =v r V 2 = v 2’ • • • ’ V m = v m 


dv l dv 2 dv m 

T m+1 = ~£f’ V m + 2 - Yf’ • • • ’ V 2m 


In terms of these variables, equations (16) can be written as a set of first-order differ- 
ential equations. Hence, we have 

~ = A fl V + s(0) (B2 


where 


2m x 1 2m x 2m 2m x 1 2m x 1 


m x m m x ml 


2m x 2m 


6 

Lm x m 


m x ml 


8 ( 9 ) = 

2m x l 


m x 1 

s(0) 

m x l 
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The solution of equation (B2) is well known and can be written as (ref. 16) 


V(0) = e Ad V(0) + e A0 / e" Ar? s(t?) dr? (B5) 

J 0 

2m x 1 2m x 2m 2m x 1 2m x 2m 2m x 2m 2m x 1 


where V(0) is an initial value vector whose evaluation from given boundary conditions is 

/v 

A 8 

discussed in appendix C and e is a matrix series given by 



i + M + aV + aV + . 

1! 2! 3! 


(B6) 


In terms of the coefficient matrix Kg, equation (B6) yields 





e 2i+1 

(2i + 1)1 


K 


oo oo 

V i 21tl „i+i 

/ , (2i + l): 9 (2i): 9 


Li=o 


i=0 


A n (0) 

A l 2 (0 ) 

II 

< 

V 

m x m 

m x m 

2m x 2m 

A 2i( 6 ) 

A 22 (0 ) 


_m x m 

m x m. 


From equations (B7) and (B8) we note that 


A n (0) = A 22 (d) = cosh Kg/ 2 ? 
A 12 (?) = Kg 1 / 2 sinh K*/ 2 ? 


(B7) 


(B8) 


(B9) 

(BIO) 
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(Bll) 


A 21 (6)=K 0 A 12 (0) 


1/2 

where Kg' is a matrix whose square is equal to Kg. Noting that and A 22 are 

even functions of 7 ) while and A,^ are odd functions of 6 , the inverse of 

e A 9 becomes 


e 




-A 6 


A n G) -a 12 (0 j! 

-A 21 (S) A 22 ( 0 )_j 


(B12) 


Substituting equations (B8) and (B12) into the identity of e^® • e~ A0 = I yields 


A^e)- Kg k\ 2 Ce)= i (Bi3) 

m x m m x m m x m m x m 


Matrix functions (B9) and (BIO) may be evaluated for each value of the independent vari- 
able from the series definitions given in equation (B7). The truncation point in each 
series is determined by the required accuracy for the convergence of the successive 
approximation calculations. Equation (B13) may be used as a measure of this required 
accuracy. However, for increasing values of 'd serious convergence difficulties may 
arise and an im practically large number of terms must be calculated. In order to avoid 
these computational difficulties, additive formulas for these matrix functions may be 
obtained from using the identity of 


A6 . A8 9 A(0«+0 2 ) 

e • e = e 


(B14) 


In terms of the submatrices A^e ), this identity yields the following equations: 


,/v /v 




An( 0 i + ® y) - A ll( 9 1^11^ 0 2^ + A 21^ 0 1^12^ 0 2^ 


(B15) 


a i 2 ( 0 i + 9 2^ - A 12^ 0 P A 11^ 0 2^ + A ll( 9 1^ A 12^ 9 2 ) 


,r<J 




(B16) 


The advantage of using equations (B7), (B15), and (B16) is that the solution of an eigen- 
value problem is not required. Alternate methods for the evaluation of these matrix 
functions using the solution of the deter minantal equation are discussed in reference 6. 
The particular integrals of equations (16) are constructed in a similar manner to 
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those listed in equations (A15) and (A16). The analogous particular integrals are 



»v 

rO 


B 3 (0) = - 

l A 12^ 

s(r?) dr] 

m x l 

m x m 

i-H 

x 

S 

B 4 tf)=j| 

r 9 A U (8> 

o 11 

s(r?) d?7 

m x 1 

m x m 

m x 1 


(B17) 


(B18) 


Assuming that s(ri) is known, these integrals can easily be evaluated. When the parti- 
tioned form of the matrices is used, solution vector (B5) for the circumferential differ- 
ential equations may be written as 



A n(0) A 12 (0) /' v(0jl |b 3 (0)\ 

a 21 Cb) a 22 (0)|v(o)| 


(B19) 


which when expanded yields equations (22) and (27). 
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APPENDIX C 


DEVELOPMENT OF INITIAL VALUE VECTORS 


Since the problems discussed in this report are two point boundary value problems, 
the initial values of both the displacements and their derivatives are usually not avail- 
able. As an example of how to obtain these initial values from given boundary conditions, 
consider the problem of figure 2 in more detail. 

From the definition of vector and the symmetry of the problem, we can imme- 
diately conclude that 

F 1 = ru|~ =0 = ° (Cl) 

l x 1 


The zero normal stress boundary condition on the surface r = b will be used in con- 
junction with equation (Cl) to evaluate the vector F 2 for this problem. Using equa- 
tion (5) gives 


u 


/>j ^ 

r =b 



(C2) 


Equation (18), taken at r = b, can be used to eliminate the vector u | ~ in equa- 
tion (C2). The resulting equation for the vector u | ~ jg can then be substituted into 
equation (21), which is also evaluated at r = b, to find the vector Fg. The results of 
these manipulations are 


f 2 = 


X n b l5 'l?Jb- n b ls! a B 1 <' E) - b 2< ? > 


X + 2G 


(C3) 


where 




2G 


,/^j . .(V , 


(X + 2G)b‘ 


^ll(b) - f2^^(b) 


(C4) 


n b " 


2G 


(X + 2G)b ‘ 


^12(b ) “ ^ 22 ^ ^ 


(C5) 


Note that by using L’Hospital's rule as r - 0 the solution vectors (18) and (21) for the 
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problem in figure 2 become 


u (0) = 0 



2 J 


(C6) 


Since for axisymmetric problems the circumferential displacements are inherently 
zero, the initial value vectors of the axial solutions will be considered next. An analo- 
gous solution to equations (22) and (27) along the z-directional lines can be written as 


|w(z? 


w(z) 


An(z ) Aj 2 ( z ) 
A 2 i(z) A 22^ z ) 


[w(0)| |B 5 (z) 

[w (0) 1 b 6 (z ) 


(C7) 


If the number of axial lines falling over the crack surface is denoted as NIC and those 
felling outside as NOC, the zero normal stress condition over the crack face and the 
symmetry condition in the crack plane result in the following equations: 


w(0) 

NIC x 1 



inside crack 


(C8) 


w(0) =0|~ =0 (C9) 

NOC x 1 outside crack 

Assuming that on the face z = L we have a uniform tensile stress of as shown in 
figure 2, the normal stress boundary condition on this plane gives 


w (L , = — ^2 (CIO) 

X + 2G X + 2G \ r/~ 

n x 1 


This vector can be suitably partitioned into vectors w Q; (L ) and w^(L). 

NIC x 1 NOC x 1 

ience of matrix manipulations, we construct the vectors as follows: 


For conven- 
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\ 

5a 



NIC x 

w(0? 


r > 
F 5 


F 5/3 

n x 1 




NOC x 


> = < 


► = ■< 


w(0) 


F 6 


F 6a 

n x 1 

J 


^ J 


NIC x 



F 6/3 
NOC x 


■\ 


y 


(Cll) 


Values of Fg^ and Fg^ are given by equations (C8) and (C9), respectively. From 




equation (C7) we can express w (L ) in a partitioned form consistent with equation (Cll) 
as 


w 


a 

NIC x 1 


w /3 

NOC x 1 


L 2lM 


L 21a2 


NIC x NIC NIC x NOC 

A 21/31 A 21/32 

z = £ L noc x NIC NOC x NOCj~ 


*5a 
NIC X 1 


°5a 
NIC X 1 




F 50 

NOC X 1 


B 5/3 

,NOC x 1 


22a 1 


22a2 


NIC x NIC NIC x NOC 


"22/31 "22/32 

NOC x NIC NOC x NOC 


/v /v/ 

' J z=L 


6a 

NIC x 1 


r 6/3 

vNOC X 1 


z =L / 

B 6a 
NIC X 1 


> + < 


B 


6/3 


vNOC x 1 


z =L 

/ 

(C12) 


Equation (C12) leads to two matrix equations involving the two unknown vectors Fg ff 
and Fg^. A solution of these equations yields 
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X + 2G 


u /3 + 


NOC x NOC NOC x 1 NOC x NOC NOC x 1 NOC x 1 


1 

X + 2G 


A 2l0!l °0a 


NOC x NOC NOC x NIC NIC x NIC NIC x 1 


- X A" 


. -1 / ~ u a 

A 21o;l ( u a + "y 


NOC x NOC NOC x NIC NIC X NIC NIC x 1 NIC X 1 


a; 1 A b B 5(3 (L) - a; 1 A c Bg a (L ) 

NOC x NOC NOC x NOC NOC x 1 NOC x NOC NOC x NIC NIC x 1 


B 6e (L) +_A_ 
X + 2G 

NOC X 1 


u a + 


NOC X NOC NOC x NIC NIC x l NIC x 1 


where 


A 21al A 22a2 )~ = 


NOC x NOC NOC x NOC NOC x NIC NIC x NIC NIC x NOC 


"( A 21/32 


A 21al A 21a2 /~ ~ 
' z =L 


NOC x NOC NOC x NOC NOC x NIC NIC x NIC NIC x NOC 


"( A 22/31 


A" ^ A 

a 21q-1 a 22q-1 


NOC x NIC NOC x NIC NOC x NIC NIC x NIC NIC x NIC 


(C13) 


(Cl 4) 


(C15) 


(C16) 


Note that the particular integrals, the applied stress vector, and the radial displace- 
ments and their derivatives are also partitioned according to their location with respect 


I 


to the crack. 

The crack opening displacement is given by 


5a 


NIC x 1 


A + 2 G 


A" ^ a 

A 21al a 0 a 

NIC x NIC NIC x 1 


,-l 

v 21al 


/v 



rv; 



r 



NIC x NIC NIC x 1 NIC x 1 


_B 5a (B) “ A 21al A 21a2 B 5/3 (B) 

NIC x 1 NIC x NIC NIC x NOC NOC x 1 


A 21al A 22a2 B 6a (B) 
NIC x NIC NIC x NIC NIC x 1 


A + 2G 


4lal 


22 a 1 








r 



NIC x NIC NIC x NIC NIC x 1 NIC x 1 


“ A 21al A 22a2 F 6/3 " A 21al A 22a2 B 6/3 (B) (C17) 

NIC X NIC NIC x NOC NOC X 1 NIC x NIC NIC X NOC NOC x 1 

where Fg^ is given by equation (C13). Although the matrix A 21 (L ) is singular, the 
partitioned matrices are not. In calculating the previous equations the indeterminate 
terms at r = 0 must be carefully considered. 
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